Transport coefficients of a mesoscopic fluid dynamics model 
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We investigate the properties of stochastic rotation dynamics, a mesoscopic model used for 
simulating fluctuating hydrodynamics. Analytical results are given for the transport coefficients. We 
discuss the most efficient way of measuring the transport properties and obtain excellent agreement 
between the theoretical and numerical calculations. 
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I. INTRODUCTION 



Investigating the hydrodynamics of complex fluids is a 
challenging task because of the large number of degrees 
of freedom involved, the interplay between the flow and 
thermodynamic properties and the existence of widely 
varying length and time scales. Examples include the 
dynamics of polymers and biopolymers in solution, poly- 
mer melts, amphiphilic systems, particulate suspensions 
and chemically reacting fluids. 

The complexity of these systems points to numerical 
approaches and it has become apparent that mesoscale 
algorithms such as dissipative particle dynamics [1], 
smoothed particle hydrodynamics [2] and lattice Boltz- 
mann algorithms [3,4] are increasingly providing useful 
numerical tools to investigate complex fluids and, in par- 
ticular, the role of hydrodynamics in their behaviour. 
These approaches are essentially novel ways of solving the 
Navier-Stokes equations designed to incorporate the most 
relevant microscopic physics of a given complex fluid. 

Recently Malevanets and Kapral described one such 
mesoscopic approach which we shall term stochastic ro- 
tation dynamics [5]. This method solves the thermohy- 
drodynamic equations of motion by following the paths of 
particles moving in discrete time but continuous space. 
The algorithm has two major advantages. Firstly, be- 
cause it is a particulate method, it is easy to couple it to 
a solute, say, polymers or colloids, which can be treated 
using a molecular dynamics algorithm [6-9]. Secondly, 
the algorithm is rather simple and therefore analytical 
work is possible which can be aimed at increasing our 
understanding of mesoscopic approaches to the hydrody- 
namics of complex fluids. 

Malevanets and Kapral mapped the stochastic rotation 
algorithm to a Boltzmann equation and showed that the 
model possesses an H-theorem [5]. Ihle and Kroll ex- 
tended the algorithm to correct for spatial correlation ef- 
fects and gave approximate results for the viscosity and 
diffusion coefficient in two dimensions [10]. Allahyarov 
and Gompper [11] used a Poiseuille flow geometry to 
mesure the viscosity in three dimensions. Stochastic rota- 
tion dynamics has been applied to colloids in solution [6] , 



dilute polymer solutions [7-9], amphiphilic systems [12], 
flow in porous media [13], binary fluid mixtures [14], clus- 
ter dynamics [15] and flow around solid objects [16-18]. 

To understand as fully as possible how stochastic ro- 
tation dynamics can be used to simulate hydrodynamic 
behaviour it is helpful to increase our understanding of 
its properties. Therefore this paper presents an inves- 
tigation of the viscosity and diffusion coefficient of the 
model in two and three-dimensions. 

The algorithm is described in section 2 following Mal- 
evanets and Kapral [5] and Ihle and Kroll [10]. We then 
show how it is possible to model a shear flow by imposing 
Lees-Edwards boundary conditions [19,20]. The ability 
to simulate shear is important as it provides a computa- 
tionally feasible way of measuring viscosities over a wide 
parameter range. It will also extend the range of appli- 
cations of the approach to study viscoelastic effects such 
as the shear-thinning of a polymer solution [21]. 

In section 3 we discuss the shear viscosity, 77, of stochas- 
tic rotation dynamics. We argue that it is useful to di- 
vide the viscosity into two contributions, a kinetic and a 
collisional term. We first provide analytical expressions 
for the kinetic and collisional parts of the viscosity. We 
then compare these to numerical results and discuss their 
relative magnitudes. 

Section 4 deals with the coupling of a test particle to 
the solvent model. A massive particle obeys Langevin 
dynamics and we obtain an expression for the friction 
coefficient describing the coupling to the solvent. We 
stress that the friction coefficient is related to the col- 
lisional part of the viscosity, not to the total viscosity, 
a distinction unimportant in a real fluid, but pertinent 
here. 

Finally a conclusion summarises the paper and sug- 
gests possibilities for future work. 



II. THE MODEL 

A. Stochastic rotation dynamics 

The solvent is modelled by a large number N , typically 
■ 10 5 , of point-like particles of mass m which move in 
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continuous space with a continuous distribution of veloc- 
ities, but discretely in time [5]. The algorithm proceeds 
in two steps. In the first of these, a free streaming step, 
the positions of the solvent particles at time t, Vi(t), are 
updated simultaneously according to 

n(t+5t) = n{t)+Vi(t)5t (1) 

where v i(t) is the velocity of a particle and St is the value 
of the discretized time step of the solvent which we take 
to be unity. 

The second part of the algorithm is the collision step. 
The system is coarse-grained into L d /a d cells of a d- 
dimensional lattice with side a. There is no restriction 
on the total number of solvent particles in each cell, (al- 
though the total number of particles in the system is 
conserved). Stochastic multi-particle collisions are per- 
formed within each individual cell, by rotating the veloc- 
ity of each particle relative to the centre of mass velocity 
v C m(t) of all the particles within that cell 

Vi (t+5t) = v cm (t) +R( Vi (t)- v cm (t) ) . (2) 

R is a rotation matrix which rotates velocities by a fixed 
angle a around an axis generated randomly for each cell 
and at each time step. 

The aim of the collision step is to transfer momen- 
tum between the particles while conserving the total mo- 
mentum and energy of each cell. Because mass, momen- 
tum and energy are conserved locally the thermohydro- 
dynamic equations of motion are captured in the con- 
tinuum limit [5]. Hence hydrodynamic interactions can 
be propagated by the solvent. Note, however, that any 
molecular details of the solvent are excluded: this allows 
the hydrodynamic interactions to be modelled with min- 
imal computational expense. 

The volume in phase space is invariant under both 
the free streaming and collision steps. Hence the system 
is described by a microcanonical distribution at equilib- 
rium. Also it has been proved that the algorithm obeys 
an H-theorcm [5]. 

A particularly useful feature of stochastic rotation dy- 
namics is the ease with which hydrodynamic interactions 
can be "turncd-off" thus replacing the hydrodynamic 
heat bath by a Brownian (random) heat bath. This is 
achieved by randomly interchanging the velocities of all 
the solvent particles after each collision step, thus re- 
laxing the constraint of momentum conservation from a 
local to a global one. Accordingly the velocity correla- 
tions which result in hydrodynamic interactions disap- 
pear from the fluid although the equilibrium Maxwell- 
Boltzmann distribution is maintained. This greatly fa- 
cilitates pinpointing the effects of the hydrodynamic in- 
teractions in a given situation [8]. 



B. Simulation parameters 

Unless otherwise stated m = 1, a = 1 and the num- 
ber of cells was L 3 — 32 3 . The initial solvent distri- 
bution was generated by assigning positions randomly 
within the system with an average number of particles 
per unit cell 7 = 5. Thus the total number of parti- 
cles was N = 163,840. The velocities were assigned from 
a uniform distribution (—v max <v a <v max ), a — x,y,z 

where v max = ^/ 3fc T ° T - The distribution relaxed rapidly 
(i~100) to the equilibrium Maxwcll-Boltzmann form cor- 
responding to the temperature T. Because of the long 
time-scales studied, any small non-zero net momentum 
in the system can eventually cause a large shift of a ve- 
locity distribution along the y axis. To overcome this we 
perform a Galilean transformation to remove any net mo- 
mentum from the system and then rescale the velocities 
to restore the correct temperature. 

C. Lees-Edwards shear boundary condition 

To achieve a steady shear flow we use Lees-Edwards 
boundary conditions [19,20], which are rather naturally 
incorporated into stochastic rotation dynamics. Our aim 
is to introduce a flow in the ^-direction, with a con- 
stant velocity gradient along y. Periodic boundary con- 
ditions are applied along x and z. Along y the boundary 
conditions are modified so that if a particle crosses the 
lower boundary at (x, y = 0, z) it returns to the upper 
boundary at the position (x = (x + ut) mod L ,y = Lz) 
with velocity (v x = v x + u, v y , v z ) where mod.L de- 
notes modulus L. Similarly if a particle leaves the up- 
per boundary at (x,y = L,z), it reappears at the lower 
boundary at (x = (x — ut) mod L ,y = 0, z) with velocity 
{v x = v x — u, v y , v z ). This algorithm gives a steady shear 
flow, and the system mimics one of infinite size. 

When an external force is applied to the system, energy 
is constantly produced. Moreover a small net momentum 
can be introduced by the Lees-Edwards boundary condi- 
tion for any system with a finite number of particles. To 
sustain a zero net momentum and temperature the ap- 
propriate shear velocity is subtracted from each particle 
velocity, a Galilean transformation is applied, the solvent 
velocities are rescaled to give a temperature T and the 
shear velocity is restored. This procedure is carried out 
after each collision step. 

D. Grid shift for periodic and shear boundary 
conditions 

Ihle and Kroll [10] pointed out that at low tempera- 
tures the transport coefficients of stochastic rotation dy- 
namics showed anomalies. This happens because at low 
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temperatures particles in a given cell can remain in that 
cell over several time steps and participate in a series of 
collisions. Accordingly the molecular chaos assumption 
breaks down. However Ihle and Kroll [10] showed that 
it is possible to improve the behaviour of the algorithm 
by placing the cubic grid in a random position at each 
collision step. 

In practice the easiest way to implement such a grid 
shift is to move all the particles with the same random 
vector with components in the interval [—a/2, a/2] before 
each collision step. After the collision step the particles 
are returned to their original positions. If the displace- 
ment takes the particles across the edge of the simulation 
box Lees-Edwards or periodic boundary conditions are 
applied as appropriate. 

We note that the grid shift accelerates momentum 
transport between cells. As a result the system relaxes 
more quickly. It must also be taken into account when 
calculating transport coefficients such as the viscosity. 



III. VISCOSITY 

One of the most important parameters characterising 
a fluid is its viscosity. In this section we investigate the 
shear viscosity associated with the stochastic rotation 
model. The viscosity is usefully divided into two contri- 
butions. The kinetic viscosity results from the momen- 
tum transferred during the free streaming step, and the 
collisional viscosity results from the momentum trans- 
ferred during the velocity rotations. 



A. Calculation of the kinetic viscosity 

To calculate the kinetic viscosity we consider a system 
undergoing shear with rate 7 = dU Qy V ^ in the x direction 
and use a kinetic theory approach. Consider first two 
dimensions. On average the velocity profile is given by 
v = (■yy, 0). In the steady state there is a frictional force 
acting on any plane in the fluid perpendicular to y which 
is given by the element of the stress tensor 



/oo />0 />oo 

dv x dy dv y v x P(v x 1U-, v y) 
-oo J — oo J—y/St 

/oo poc p—y/St 

dv x dy dv y v x P(v x - jy, v y 

-00 J J —oo 



"xy — st 



>1- 



d u x (y) 
dy 



(3) 



where 77 is the shear viscosity. This is the experimental 
definition of the viscosity [20]. Our aim is to calculate 
the stress tensor during the streaming step: a xy = —(flux 
of a;- momentum crossing a plane of constant y) . 

For simplicity we take y = y = 0. Consider a parti- 
cle at position (x,y) with velocity (v x , v y ). Particles will 
only cross the plane during a time step if they are moving 
towards the plane and have a velocity such that \v y 5t\ is 
greater than the distance to the plane. Therefore, the 
stress tensor can be written 



where P(v x ,v y ) is the velocity probability distribution of 
particles in the rest frame of the fluid and p = is the 
mass density. 7 denotes the average number of particles 
per cell and d = 2,3 is the dimensionality. By making 
the change of variable v' x — v x — jy, and changing the 
order of integration this reduces to 



xv — 



7/9 St 



(v y )-p{v x v y ) 



(4) 



where the averaging is over the probability distribution 
P. The first and second terms result from the shear pro- 
file and from induced correlations between v x and v y re- 
spectively. In the steady state the velocity distribution 
P deviates from the Maxwell Boltzmann form because of 
these correlations. 

To find the behaviour of (v x v y ) we calculate separately 
the effect of the streaming and collision operations. Con- 
sider the velocity distribution at the plane y = yo. At 
time t, at yo, particles with positive velocity v y come from 
2/o — VySt due to the streaming step, and correspondingly 
have a smaller x velocity. Conversly particles from y > y 
tend to have a higher x velocity. This means that the ve- 
locity probability distribution is sheared by the stream- 
ing. That is P a f ter (v x ,v y ) = P be f° re (v x +jv y St,v y ). 
Averaging v x v y over this new distribution and changing 
variables gives 



/ \after 
{VxV y ) 



/OO />00 
dv x / dv y v x v y P(v x + A /v y St,v y ) 
- 00 J — 00 

{v x Vy) - j5t (v 2 y ) . (5) 



Thus the streaming operation changes the value of (v x v y ) 
by — jSt (v y ), and tends to make v x and v y increasingly 
anticorrelated. 

Next we consider the effect of the collision step on the 
velocity distribution. Collisions tend to reduce correla- 
tions. To see this we consider a cell of n particles and 
the two dimensional collision equation 

, x , ( cosa ±sina \ . , . . 
V ( t + S V={ T sina cosa J K*) " »cm) + ^ 

where the sign corresponds to the direction of the ro- 
tation axis. The centre of mass velocity v cm can be 
divided into two contributions, v cm = (v + v) /n, the 
first from the test particle of velocity v and the sec- 
ond v = Y^i=i v i fr° m the n — 1 particles sharing the 
same cell. We assume molecular chaos, that is the ve- 
locities of different particles are uncorrelated. Then, 
(v x v y ) = (v x v y ) = and (v x v y ) = (n — 1) {v x v y ) and 
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(V X (t + St)Vy(t + St)) 

n — 1 



(1 — cos 2a) 



(v x (t)Vy(t)) . 



where we have averaged over the direction of the rotation 
axis. In the simulations n is not constant. Therefore, we 
must consider density fluctuations. The probability of 
n particles being in a given cell is given by the Poisson 
distribution P (n) = *—3— . The probability of a given 
particle being in a cell with a total of n particles is " p ^") ) 
so taking an average over this distribution 



n=l ^ 



"7 



(1 — cos 2a) 



(V x {t)Vy(t)) 



f(a, 7 ) (v x (t)v y (t)) 



(6) 



During one time step we find that (v x v y ) is first 
reduced by the streaming operation and then multi- 
plied by a constant factor in the collision operation. 
In the steady state it will therefore oscillate between 
two values and we have the self-consistency condition 
{{v x v y ) - jSt (v 2 y )) f = (v x v y ). Thus 



Substitution into (4) gives 

<T xy = fTfSt (V y ) ^- + — - 



(7) 



(8) 



Using equipartition of energy, eqn. (6) and the definition 
of viscosity (3) yields 



2D _ -yk B T5t 

Vkin 2 



(9) 



(7 - 1 + e-T)(l - cos 2a) 
Applying the same proccedure in three dimensions gives 



4Z = 

-yk B T5t 



57 



1 



(7 - 1 + e-T)(4 - 2cosa - 2cos2a) 2 



(10) 



B. Calculation of the collisional viscosity 

The collisional viscosity Vcol can also be calculated us- 
ing a kinetic theory approach. Consider a cubic cell with 
side a. As in the case of the kinetic viscosity a shear 
rate 7 = dW Q^ is applied along the a;-axis. The system 
is divided into two subcells by a plane at y = yo where 
< yo < «■ Let the subcell at yo < y < a contain m parti- 
cles with an average collective flow velocity u\ and that 



at < y < yo contain n 2 particles with an average velocity 
u 2 . We have u lx - u 2x = Kl j" 2 " 2 {u lx - v cmx ) where v cmx 
denotes the x-component of the centre of mass velocity 
of the cell. As the average distance between the subcell 
velocities is 5y = § the shear rate is 

du x _ uix - u 2x 
dy 



7 = 



5y 
2n 



m) 



(11) 



where n = n\ + n 2 is the number of particles in the cell. 

Our aim is to calculate the momentum crossing the 
plane at yo- We consider the n\ particles at yo < y < a. 
The momentum transfer between the two subcells 



^Pix(t + St) - ^Pl. 




(12) 



where i runs over the ni-particles, is calculated using 
the collision operation (2). Averaging over an isotropic 
distribution of the rotation axis gives 



,d-l 



St 



-ni (1 - cos a) (ui x - v c 



c) 



(13) 



where m is the mass of the solvent particle. Using eqns. 
(11), (13) and the definition of the viscosity (3) we obtain 



Vcol 



mn\ (n — ni) 
dnSt 



,d-2 



(1 — cos a) 



As n\ + n 2 is in general small we must again consider 
fluctuations in the particle density. The numbers in the 
subcells, n\ and ni are binomially distributed and aver- 
aging over them gives 



Vcol 



m (1 — cos a) 



,d-2 



dSt 



Here we used ^ni(n — n\)P{n\)un = n2 P — npq 
where p = 1 — ^ and q = Using the Poisson distri- 
bution to average over n 
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Vcol = 



m (1 — cos a) 
a d - 2 dSt 

m (1 — cos a) 



e -7 7 » 



n=2 



,d~2 



dSt 



(7- l+e~ 7 ) 



where the summation runs from 2 to 00 because there is 
no momentum transfer unless ni>l and n 2 >l. Finally, 
averaging over all planes < y < a gives 

m (1 — cos a) 



Vcol = 



(7 - 1 + e~ 7 ) 



(14) 



6a d ~ 2 dSt 

In the large density limit the effect of the density fluc- 
tuations vanish and the theory agrees with the expression 
obtained by Ihle and Kroll [10] for the two dimensional 
case (d = 2). Also for case a = ^ and d = 2 the re- 
sult presented here is identical to that given in [5] by 
Malevanets and Kapral. When d — 3 we find that our 
expression differs from that given in [5] for a = ^. 
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C. Measuring the shear viscosity numerically 

The computationally most efficient approach to mea- 
sure the shear viscosity directly exploits the experimental 
definition of viscosity (3). A shear is applied along x set- 
ting up a velocity gradient along y. The stress tensor a xy , 
equal to the flux of x-momentum across a plane perpen- 
dicular to the y axis can easily be measured numerically. 

Again we consider the momentum transfer in two 
parts, that due to the streaming step and that due to the 
collision step. A subtlety that needs to be taken into ac- 
count occurs because of the finite size of the grid imposed 
in the collision step. This means that the streaming and 
collisional contributions to the momentum transfer de- 
pend on the position of the plane chosen for the mea- 
surement. For example, if the plane lies at the edge of 
a cell, there is no momentum transfer during the colli- 
sion step whereas if the plane is taken at the centre of 
a cell the momentum transfer in the collision is max- 
imised. The streaming contribution also varies, in such a 
way that the total momentum transfer in one simulation 
time step is independent of the position of the measuring 
plane. This is illustrated in figure 1. To separate out 
contributions to the kinetic and collisional viscosities it 
is necessary to average over all plane positions within a 
cell: this is another advantage of imposing the grid shift 
which performs the averaging over the cell automatically. 

To measure the shear viscosity we used a simulation 
of size 32 x 32 x 32. Four thousand time steps were re- 
quired to reach the steady state and the stress tensor 
was averaged over ten thousand time steps. It is also 
necessary to average over many measurement planes to 
minimise thermal noise. We choose 25 planes from the 
32 cubed box. Planes near the edge of the simulation 
box are omitted to avoid edge effects resulting from the 
Lees-Edwards boundary conditions. We also note that 
the pressure, P = —<y yyi can be measured in the same 
way, and takes the expected value for an ideal gas. 

To check that the algorithm corresponds to a Newto- 
nian fluid we plot in figure 2 the dependence of the shear 
viscosity on the shear rate. The viscosity is independent 
of the shear rate, as expected, except at very low and 
very high shear rates. The deviations at low shear are 
due to statistical errors because the shear represents only 
a small perturbation to the Maxwcllian velocity distribu- 
tion. Those at high shear result from finite-size effects. 
They occur when the distance moved by the walls in each 
time step approaches the size of the system. The value 
of 7 chosen to measure the shear viscosity was 0.0625 at 
k B T = 0.8. 

The shear viscosity t]gk can also be measured using an 
equilibrium approach based on the Green-Kubo formula 
[6,10,22] 

V gk= hm / dS(a xy (S)<T xy (0)). (15) 



For discrete time steps this can be approximated by 

Vgk = 

VSt ( 1 °° \ 

■j^j ( 2 <^w(°)^(°)> + E (^y(^t)a xy (0))j (16) 

where St is the time step of the solvent dynamics. 

Further details of the equilibrium approach is de- 
scribed in [10]. However thermal fluctuations mean that 
very long runs are needed to obtain statistically mean- 
ingful results - using shear is faster by a factor ~ 30. 
Therefore this approach was used only as a check on the 
validity of the non-equilibrium method. 

D. Results for the viscosity 

Consider first the kinetic viscosity. Figures 3 and 4 
present the dependence of the three dimensional kinetic 
viscosity on the rotation angle a at two different tem- 
peratures. Agreement between the theoretical result eqn 
(10) and the simulations is excellent. As a goes to zero 
the kinetic viscosity diverges. This corresponds to parti- 
cles moving with an infinite mean free path, unimpeded 
by collisions. From eqn (6) and (7), as a — ► the cor- 
relation (v x v y ) — > — oo. In two dimensions the kinetic 
viscosity is symmetrical about a = \ and hence also 
diverges at a = tt. This differs from the three dimen- 
sional case where the viscosity reaches a local maximum 
at a = tt. 

A small deviation between the theory and the numer- 
ical results is observed at low temperatures. We believe 
that this is due to one of two effects. At such low temper- 
atures, even with a gridshift, the assumption of molecular 
chaos is an approximation. Also at these temperatures 
the equilibrium distribution no longer approximates to 
a Maxwell-Boltzmann distribution due to the effect of 
the shear. Figure 5 compares the density dependence of 
the kinetic viscosity with theory, again showing excellent 
agreement. 

We now turn to the collisional viscosity. This is as pre- 
dicted, independent of temperature (eqn (14)). Figures 6 
and 7 show its dependence on rotation angle a and par- 
ticle density respectively. Agreement with the theory is 
again excellent except at extremely low temperatures. 

The kinetic viscosity is proportional to the tempera- 
ture and we recall that the collisional viscosity does not 
depend on the temperature. Hence the former dominates 
at high temperatures, and the latter at low temperatures. 
This has important indications for the balance of viscous 
and diffusive transport as discussed below. 

Finally we check for any finite size effects in the sim- 
ulations. Figure 8 shows that for the systems studied 
there is only a very small effect on the viscosity at the 
smallest system size. 
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IV. COUPLING TO A SOLUTE PARTICLE 



B. Numerical results for the friction coefficient 



To place solute particles such as polymers or colloids 
in the stochastic rotation dynamics solvent it is impor- 
tant to understand the coupling between them and the 
solvent. Therefore we calculate here the friction coeffi- 
cient acting on a particle of mass M and velocity v in 
a Brownian solvent, and then show that the fluctuation- 
dissipation theorem holds as expected. 



A. Friction coefficient 

The discretized Langevin equation for a monomer is 
r v(t+5t)-v(t) 



M- 



5t 



-£v(t) + v(t) 



(17) 



where £ denotes the friction and St is the time step. 77(f) 
is stochastic noise with time average (77) =0. The evo- 
lution of the velocity is given by eqn (2) which may be 
rewritten 

v(t + St) - v{t) = {R-I) (v(t) - v cm ) (18) 

where R is the usual rotation matrix and I is the identity 
matrix. 

We consider a particle of mass M in a cell with n sol- 
vent particles of mass to. The centre of mass velocity of 
the cell is given by 



Mv(t) + J27=i mu i 

M + mn 



(19) 



Averaging eqn (18) over an isotropic distribution of the 
rotation axis and over Ui, assuming (m) = 0, gives 

2 711 T) 

(v(t+St) v(t)) = — (1 - cos a) <„(*)> . (20) 

In time the number of particles within a given cell, n, 
changes. The probability is given by the Poisson distri- 
bution, P(n) = e , where 7 is the average number per 
cell. By comparing (20) with the time averaged Langevin 
equation (17), we identify the friction coefficient as 

,sn e~ 7 7™ mn 2M . . . , 

£ 3 = E J U , — TT. (! " cos «) • ( 21 ) 

In two dimensions 



n\ M + mn 3 St 



e D = E 



n=0 



e 7 7 n mn M 
n\ M + mn St 



(1-cosa). (22) 



The friction coefficient is measured numerically using 
eqn (17). Multiplying (17) by v(t) and taking a time 
average 

(v(t+St).v(t))-(v 2 (t)) = -^(v\t)) (23) 

where we use (77(f) •«(£)) = (f](t))(v(t)) = 0. Therefore, 
using the cquipartition theorem 



M ( M 



(v(t+St)-v(t)} 



(24) 



Theory and numerical results are in good agreement for 
different M as shown in figure 9. 



C. Fluctuation-dissipation theorem 

Next we verify fluctuation-dissipation theorem. Sub- 
stituting (18) into (17) and rearranging gives 

77(f) =^(R - I) (v(t) - v cm ) + £„(t). (25) 

To calculate the equal time correlation function we must 
assume equipartition of energy and molecular chaos, 
therefore {v a Uifj) = (v a ) (uip) = and {ui a Ujp) — 
SijSapkBT /m. Using eqn (19) and averaging over an 
isotropic distribution of the rotation axis and the solvent 
velocity distribution we find 

(77(f). 77(f)) = M ^ BT ( 1 - ™ { }- C ° S " } ) (26) 



St 



d(M + mn) 



for dimension d. The first term in the bracket is that 
obtained from the fluctuation-dissipation theorem. The 
second is a small correction, which goes to zero in the 
limit of large M. 



D. Diffusion coefficient 

The diffusion coefficient can be calculated using the 
definition 



D = lim 

t—*oo 



(r(t)-r(O))' 

Mt 



1 1 00 

- (v (0) • v (0)) St+-J2(v (nSt) ■ v (0)) St. (27) 



2d 



n=l 



We note that the resulting equation (27) is the same as 
the Green-Kubo formula, which is derived from linear 
response theory [22]. 

For a Brownian particle [23] the velocity-velocity au- 
tocorrelation is easily calculated using eqn (23) 



G 



(v(t+nSt)-v(t)) = (l- ^J(v 2 (t)) . (28) 

Substitution into (27) gives, 

k B T ( 8tf\ 

This is the normal Einstein relation with an error term 
which becomes vanishingly small in the large M limit. 

V. CONCLUSION 

We have provided analytic expressions for the shear 
viscosity of a two- and three-dimensional mesoscale 
model of hydrodynamics, stochastic rotation dynamics. 
The derivation assumes molecular chaos. The viscosity 
has two contribution r)ki n and i] co i from the streaming 
and collision steps of the algorithm respectively. The 
former dominates at high temperatures and the latter at 
low temperatures. 

We also obtained an expression for the viscous drag £ 
acting on a test particle. Note that £ and r] co i given by 
eqns (21) and (14) have the same dependence on the ro- 
tation angle a and are both independent of temperature. 
This is to be expected as they both result from momen- 
tum transfer in the collision step of the algorithm. It is 
useful to realise that in the stochastic rotation method 
for a fluid with a given total viscosity the coupling to a 
test particle will depend on temperature and will be very 
weak at high temperatures where r\un r\ C ol- 

The viscosity in three dimensions was measured nu- 
merically by using Lees-Edwards boundary conditions to 
apply a shear to the system. The analytic and numeri- 
cal results were found to be in excellent agreement thus 
providing a check on both calculations. It is important 
to have a robust way of measuring transport coefficients 
for this model as it can easily be coupled to a solute such 
as polymers or colloids where the non-Newtonian flow 
properties of the resultant solution are of interest. 
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FIG. 1. Dependence of the shear viscosity r\ on the position of the measurement plane relative to a cell for feT = 0.8, 
7 = 5.0 and a = \. 
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FIG. 2. Dependence of the shear viscosity ij on shear rate 7 for ksT = 0.8, 7 = 5.0 and a = ^. 
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FIG. 3. Dependence of the kinetic viscosity r\uin on rotation angle a for UbT = 0.8 and 7 = 5.0. 
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FIG. 4. Dependence of the kinetic viscosity r/tin on rotation angle a for feT = 0.04 and 7 = 5.0. 
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FIG. 5. Dependence of the kinetic viscosity r] k i n on mean particle number 7 for fcsT = 0.8 and a = f. 
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FIG. 6. Dependence of the collisional viscosity r\ co i on rotation angle a for 7 = 5.0. 
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FIG. 7. Dependence of the collisional viscosity r; coi on mean particle number 7 for fesT = 0.8 and 
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FIG. 8. Effect of finite lattice size on the measured viscosity for fcsT = 0.8, 7 = 5 and a — f 
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FIG. 9. Dependence of the friction coefficient £ on rotation angle a for different particle masses M for fcsT = 0.8 and 7 = 5.0. 



12 



